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STATISTICS  OF  COMPLEX  CROSS  SPECTRUM  ESTIMATE  FOR 
SINUSOIDAL  SIGNALS  AND  ARBITRARY  NOISE  SPECTRA 

INTRODUCTION 

The  analysis  of  the  stability  of  the  estimate  of  the  complex 
cross  spectrum  usually  proceeds  on  the  basis  that  the  two  input 
processes  have  a  slowly  varying  cross  spectrum  relative  to  the 
spectral  window  employed.  See,  for  example,  [1;  (4)  and  sequel]. 
Here,  we  will  eliminate  that  restriction  and  allow  real  input 
signals  with  arbitrarily  narrow  width,  namely  sinusoids,  and 
allow  real  additive  input  noises  with  arbitrary  spectra. 

On  the  other  hand,  we  will  restrict  consideration  to  the 
special  case  where  the  two  input  noise  processes  are  zero  mean 
Gaussian  and  are  statistically  independent  of  each  other. 
Furthermore,  the  temporal  weightings  applied  will  be  presumed 
nonoverlapping  in  time,  thereby  leading  to  (approximately) 
independent  spectral  estimates  for  each  time  segment. 

We  will  derive  the  exact  joint  characteristic  function  of  the 
real  and  imaginary  parts  R  and  Q,  of  the  cross  spectrum  estimate 
G  »  R  -i-  iQ,  with  both  signal  and  noise  present.  This  result 
enables  determination  of  the  high-order  joint  cumulants  of  random 
variables  R  and  Q,  for  arbitrary  signal-to-noise  ratios. 

For  the  noise-only  case,  the  corresponding  joint  probability 
density  function  of  random  variables  R  and  Q  will  be  derived  in 
closed  form.  It  can  then  be  used  to  derive  various  moments  of 
the  magnitude  of  estimate  G,  such  as  the  average  magnitude  |g|. 
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When  signal  is  present,  the  joint  probability  density 
function  of  R  and  Q  cannot  be  found  in  closed  form.  Instead,  a 
two-dimensional  fast  Fourier  transform  is  utilized,  followed  by 
numerical  integration  to  find  the  moments  of  interest.  Compari¬ 
son  of  these  accurate  results  with  a  Gaussian  approximation 
affords  quantitative  'verification  of  the  Gaussian  approximation 
when  the  number  of  pieces,  K,  used  in  the  finite  average  for 
cross  spectrum  estimate  G  exceeds  10  approximately. 

A  deflec  ion  criterion  of  cross  spectrum  magnitude  estimate 
|6|  is  defined  and  evaluated,  both  numerically  and  by  use  of  the 
same  Gaussian  approximation  for  che  joint  probability  density 
function  of  R  and  Q.  Finally,  the  same  statistics  are  evaluated 
for  an  auto  spectrum  estimate  obtained  from  the  two  input 
processes  optimally  scaled  prior  to  addition. 


2 


TR  10709 


PROBLEM  DEFINITION 

This  study  is  a  follow-on  to  an  earlier  report  [1]/  where  the 
overlapped  fast  Fourier  transform  processing  method  of  weighted 
data  segments  for  the  purpose  of  estimation  of  the  cross  spectrum 
was  well  documented.  Familiarity  with  that  material  and  results 
is  presumed  in  the  following  development. 

A  common  sinusoidal  signal  s(t)  at  frequency  f^  is  present  in 
both  input  channels;  that  is, 

s(t)  -  A  cos(2iif^t  +  f)  ,  (1) 

where  f  is  a  random  variable  uniformly  distributed  over  2ii.  The 
input  noise  processes  in  the  two  channels  are  x{t)  and  y(t), 
respectively.  The  k-th  time  weighting  function  and  its  spectral 
window  (Fourier  transform)  are  given  by 

w^(t)  -  w(t  -  Tj^)  for  1  i  k  1  K  ,  (2) 

Wk(v)  s  J  dt  exp(-i2nvt)  Wj^(t)  =  exp( -i2rvTj^)  W(v)  ,  (3) 

where  K  is  the  total  number  of  pieces  used  in  the  cross  spectnun 
estimate  G,  and  W(v)  is  the  window  (Fourier  transform) 
corresponding  to  fundamental  time  weighting  w(t).  Time  delays 
|Tj^)  are  taken  widely  enough  spaced  that  individual  temporal 
weightings  {Wj^(t) }  do  not  overlap  on  the  time  axis.  (Integrals 
without  limits  are  over  the  range  of  nonzero  integrand.) 

The  k-th  voltage  density  estimate,  at  analysis  frequency  f, 
of  the  signal  component  is,  for  both  channels. 
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Sj^  ■  J  dt  exp(-i2nft)  Wj^(t)  A  cos(2nf^t  +  ♦)  *  (4) 

a  I  J  dt  w^(t)  exp(-i2n(f-£^)t)  -  |  W^(f-fQ)  =  a,^  +  / 

where  random  variables  and  are  real.  It  is  presumed  that 
analysis  frequency  f  is  close  to  the  actual  input  signal  center 
frequency  but  they  need  not  be  equal. 

The  k-th  voltage  density  estimate,  at  analysis  frequency  f, 
of  the  noise  component  x(t)  is 

Xr  *  J  dt  exp(-i2nft)  Wj^(t)  x(t)  =  aj^  +  ibj^  ,  (5) 


for  1  ^  k  1  K.  Random  variables  aj^  and  bj^  are  zero  mean 
Gaussian,  since  the  two  input  processes  x(t)  and  y(t)  are  zero 
mean  Gaussian  processes.  The  corresponding  quantities  for  the 
other  channel  are 


Yr  «  J  dt  exp(-i2iift)  Wj^(t)  y(t)  *  Cj^  +  • 


(6) 


The  complex  cross  spectrum  estimate  at  analysis  frequency  f 
is  therefore  given  by 


G  ■  ^  ZZ  (Sjj  +  *k)<\  +  = 


1 

K 


k-1 


(a,,  +  iPv  +  «!,  +  ib,  )  (a.  -  ip^  +  c.  -  id.  ) 


1 

K 


K 

m  (Rk  +  iQk)  s  R  +  iQ  / 
k-1  ^ 


where  real  and  imaginary  parts 


(7) 


4 


TR  10709 


Qk  -  (P|j  +  b^)(Ok  +  C|t)  -  (Ojt  *k>^^k  '*■  **k)  * 

It  is  important  to  notice  from  (4)  and  (3)  that 

l®kl^  ■  l“k  ■  “k  *  ^k  *  ?  |W(f-fo)|^  ■  Y  ,  (9) 

where  the  latter  quantity  y  is  a  constant,  not  a  random  variable, 
and  that,  furthermore,  y  is  independent  of  k,  the  segment  number. 


STATISTICS  OF  aj^  AND  bj^ 

From  ( 5 ) ,  we  observe  that  ensemble  average 

l*kJ^  “  l®k  ^^kl^  “  ®k  *»k  ' 

-  JJ  dt  du  exp(-i2nf  (t-u) )  Wj^(t)  w*(u)  Rj^(t-u)  » 

-  J  dv  G^(v)  |Wj^(f-v)|^  “  J  dv  G^(v)  lW(f-v)|2  ,  (10) 

where  Rjj(t)  and  G^(v)  are,  respectively,  the  covariance  and 
spectrum  of  random  process  x(t),  and  we  used  (3).  Again,  notice 
that  this  average  is  independent  of  k.  Also,  the  result  in  (10) 

holds  regardless  of  the  relative  widths  and  variations  of  window 

2 

|W(v)|  and  spectrum  G^^Cv). 

At  the  same  time,  from  (5)  and  (3), 
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*k  “  ^“k  “  *k  ■  **k  *k  ^k  “ 

-  J  dv  G^(v)  Wj^(f-v)  Wj^(£+v)  - 

-  exp(-i4nfTj^)  J  dv  G^(v)  W(f-v)  W(f+v)  *  0  ,  (11) 

if  analysis  frequency  f  is  not  near  zero  frequency.  Combining 
(10)  and  (11)/  we  find  properties 

*k  “  ‘^k  “  M  |W(f-v)|2  s  -  0  /  (12) 

which  are  independent  of  k.  Thus,  Gaussian  random  variables  aj^ 
and  are  statistically  independent  of  each  other. 


STATISTICS  OF  Cj^  AND  dj^ 

In  an  entirely  similar  fashion,  but  working  instead  from  (6), 
^  ^  ® -  0 ' 

°k  "  ’*k  ■  M  ■  "y  '  <“> 

all  quantities  being  independent  of  k.  Furthermore,  the  time 
delays  (Tj^)  in  (2)  are  widely  enough  separated  that  all  the 
random  variables  for  weighting  k  are  independent  of  all  those  for 
weighting  m,  when  k  m. 
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JOINT  CHARACTERISTIC  FUNCTION  OF  R  AND  Q 

The  joint  characteristic  function  of  the  k-th  pair  of  random 
variables  and  Q|^  in  (7)  and  (8)  is  given  by  ensemble  average 

fj5(2/y)  ■  exp(zRj^  +  yQj^)  -  exp  z ( )  ( Oj^+Cj^ )  + 


where  we  have  let  variables 

z  ■  i5  and  y  «  i^l,  t  and  1]  real  ,  (15) 

for  shorthand  purposes.  At  the  same  time,  from  (12)  and  (13  , 
the  joint  probability  density  function  of  aj^,  bj^,  Cj^,  dj^  is, 
for  all  k,  given  by 


-1 


-1 


p(a,b,c,d)  »  (2na^]  (2«<yy)  exp 


+  b" 


+  d" 


2  a" 


2  a. 


(16) 


When  we  employ  this  result  for  p  in  the  average  required  by 
(14),  holding  random  variables  and  fixed  for  now,  the 
resulting  four-fold  integral  can  be  evaluated  by  first  evaluating 
the  double  integral  on  a  and  b,  followed  by  the  double  integral 
on  c  and  d,  by  means  of  the  following  result: 


for  Oj. 


JJ  dx  dy  exp|-  ^ouc^  -  ^Py^  +  yxy  +  +  vyj  » 

2  2 

Bu  +  gy  +  2Yt/v 
2(ap-Y^) 

2 

>  0,  p^  >  0,  cij.  Pj.  ^  Yr*  result,  after  much 


(17) 
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manipulation,  is  the  joint  characteristic  function  of  Rj^  and  Qj^, 
conditioned  on  given  fixed  values  of  random  variables  and 
namely 


£|t(2/y)  - 


,  2  2.  2  2. 

1  -  y  ) 


exp 


2  +  (2^  +  y^) (a^  +  Oy)/2 


2  2  2  2 
1  -  o‘  a‘(z^  +  y^) 


(18) 


where  we  ust  d  relation  ( 9 ) . 

But,  since  (18)  contains  no  random  variables,  it  is  actually 
the  unconditional  joint  characteristic  function  of  Rj^  and  Qj^. 
Also,  since  the  constant  y  is  independent  of  k,  we  see  that 
fj^(z,y)  is  independent  of  k.  This  leads  to  the  desired  result, 
namely  the  joint  characteristic  function  of  summation  variables  R 
and  Q  in  ( 7 ) ,  as 


fRQ(2/y)  -  exp(zR  +  yQ)  -  exp(|  HZ  *^k  K  ^  ®k) 

}c*l  Ic*  1 


1  -  ol  +  y^)/K^]  ejtp 


xl 

-  f, 

Ik'kJ 

-K 

f  2  +  (2 

1  2  2 , „2  .  „2 , /„2 
+  y  )/K 


(19) 


This  is  an  exact  result  for  the  joint  characteristic  function  of 
R  and  Q,  under  the  conditions  cited  above,  such  as  disjoint 
temporal  weightings  {Wj^(t)},  common  sinusoidal  signal  s(t),  and 
independent  input  noise  processes  x(t)  and  y(t).  The  complex 
cross  spectrum  estimate  is  given  by  (7)  as  G  ■  R  +  iQ. 
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JOINT  CUMULANTS  OF  R  AND 

When  we  expand  In 
y  according  to 


Q 

(z,y)  in  (19)  in  a  power  series  in  z  and 


In 


m»0  n«0 


^mn 


m  n 
m!  n! 


'■oo 


(20) 


then,  X|nn  is  the  m,n-th  joint  cumulant  of  R  and  Q.  There  follows 


XlO  •  Y 


^01  *  °  ^  ^mn  “  °  for  m  +  n  -  3,5,7,... 


2  0^ 

^20  “  ^02  “  I  ^  ' 


ai 


0  , 


ii  4  4 

X22 - ^  (1  ♦  *  2Ry)  , 


^40  “  ^04  “  ^  ^22  '  X31  *  Xi3  “  0  / 


(21) 


where  we  have  defined,  with  the  help  of  (9)  and  (12), 


*  2al  J  d.  Gjj(v)  |W(£-v)|^ 


i  A^  |W(£-f„)|^ 


y  2ffy  J  dv  G  (V)  |W(£-v)| 


2  * 


(22) 


These  latter  quantities  are  measures  of  the  signal-to-noise  power 
ratios  at  the  outputs  of  window  |W(v)|  in  (3). 
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The  results  in  (21)  indicate  that  as  the  number  of  pieces 
K  w,  the  two  random  variables  R  and  Q  tend  to  joint  Gaussian. 
For  example,  we  find 


X40 

X20 


3  1  ^  2R^  ^  2R^ 

(1  +  R*  **■ 


-»  0  as  K  -»  ®  . 


(23) 


Since  the  joint  third-order  moments  are  all  zero,  this  result 
indicates  a  rather  rapid  approach  to  the  Gaussian  approximation. 
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GAUSSIAN  APPROXIMATION 


The  quantities  x^q  ^01  means  of  R  and  Q, 

respectively.  Also,  X20  ^02  respective  variances, 

2 

which  are  equal  and  will  be  denoted  by  a  ;  that  xs. 


2 

o 


2  a" 


(1 


"y' 


(24) 


Since  covariance  Xj^^  is  zero,  the  Gaussian  approximation  to  the 
two-dimensional  probability  density  function  of  R  and  Q  is  given 
by 

2  2 

p  (u,v)  s  ■  ^  2  ®*p(“  ~  '^^2 

^  2nc'^  '  2o^  ^ 


MEAN  MAGNITUDE  OF  G 

With  this  approximation  at  hand,  we  can  now  evaluate  some 
moments  of  G  that  are  not  available  directly  from  joint 
characteristic  function  (19)  or  joint  cumulants  (21).  In 
particular,  we  are  interested  in  the  mean  magnitude  of  complex 
cross  spectnun  estimate  G  =  R  +  iQ;  namely, 

Pj  *  |Gl  «  (R^  +  Q^)^  =  JJ  du  dv  (u^  +  v^)**  Pj^q(u,v)  «  (26) 

-  JJ  du  dv  (u^  +  v^)**  Pg(u,v)  s  »  (27) 
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du  dv 


;u2  +  v^,*' 


2  .  „2, 


2  .  2, 


J  ^  «p(-  ^cH] 


h  h 

(5)  a  exp(-V)  jFj(1.5;l;V)  -  [|)  o  ( - . 5 ; 1 ; -V)  ,  C 


where  we  used  [2;  6.631  1],  [3;  13.1.27],  and  defined 


R..  R 

_ X  y _ 

1  ♦  "x  *  "y 


(29) 


The  function  ^Fj^  is  the  confluent  hypergeometric  function 
[3;  13.1.2  and  13.1.10]. 

It  must  be  repeated  that  (28)  is  an  approximation  for  the 

desired  average  (j^,  because  the  Gaussian  density  function 

p  (u,v),  employed  in  (27)  and  the  sequel,  is  itself  an 
y 

approximation  to  the  true  (unknown)  probability  density  function 
Pj^q(u,v)  of  R  and  Q.  The  result  of  the  Gaussian  approximation  in 
(27)  and  (28)  has  been  denoted  by 

Upon  use  of  (24),  the  mean  value  approximation,  (28), 

takes  the  form 


‘^x  °y  (t)  **x  ^y>^  iF^(-.5;l;-V)  , 


(30) 
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where  the  various  parameters  are  defined  in  (12),  (13),  (22),  and 
(29).  If  the  signal  is  absent  at  the  input,  then  A  =  0,  Rjj  * 

Ry  “  0/  V  »  0,  giving 


>-lg(A-0) 


(31) 


-k 

which  decays  to  0  as  K  ^  for  large  K. 

On  the  other  hand,  for  A  >  0,  suppose  that  K  is  large  enough 
that  parameter  V  in  (29)  is  large  compared  with  1.  Then,  we  have 
the  asymptotic  result  [4;  A. 1.16b] 


''ig  ■  *  w)  "  ?  *  4?)  AS  V  ,  (32) 

where  we  used  (9).  That  is,  the  mean  magnitude  approaches 

the  signal-only  output  y,  with  an  additive  tenn  that  decays  as 
-1  -k 

K  ,  not  K  ^  as  in  (31).  These  results  are  expected  to  be  most 
accurate  for  large  K,  where  the  Gaussian  approximation  is  best. 


MEAN  SQUARE  MAGNITUDE  OF  G 

The  mean  square  magnitude  of  cross  spectrum  estimate  G  is 


_  _  _  2  2 

/i2  *  |G|2  =  r2  +  q2  =  r  +  +  Q  +  +  20^  ,  (33) 

where  we  used  (21)  and  (24).  Upon  additional  use  of  (22),  this 
develops  into 


.22 

"2  ■  ^  "x  "y 


R  R  + 
X  y 


^  ^  ^x  ^ 
K 


)• 


(34) 
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This  result  is  exact,  having  been  developed  directly  from  the 
exact  joint  characteristic  function  (19)  of  R  and  Q;  hence,  there 
is  no  need  to  add  subscript  g  to  /J2*  However,  it  can  be  noted 
that  use  of  the  Gaussian  probability  density  approximation  (25) 
yields  exactly  the  same  result  ( 34 ) . 


DEFLECTION  OF  MAGNITUDE  |g| 

The  deflection  of  the  magnitude  of  the  complex  cross  spectrum 
estimate  |G|  is  defined  here  as 


^1  “ 

(p2(A=0)  -  /i2(A=0)]**  ‘ 


(35) 


However,  since  only  the  approximate  result  is  available,  we 
also  define  the  Gauss ian-approximation  deflection  as 


-  Pi„(A=0) 


la _ ZISL 


0)  -  ^tg(A»0 


(36) 


Substitution  of  (30),  (31),  and  (34)  (with  A  =  0)  yields 

<*9  ■  (4^)  W  "y*'’  l''l(-5il;-V)  -  1]  -  (37) 

where  V  is  given  by  (29). 

If  the  number  of  pieces  K  is  so  large  that  V  >>  1,  use  of  the 
asymptotic  behavior  of  ^Fj^  (4;  A.  1.16b]  yields 
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-  (4^  *  "x  ■'y)’’  -  (4^)’’  -  2-14  (K  R,  Ryl**  -  1- 


91  as  K  ->  »  . 

(38) 


L 

Thus,  the  Gaussian  deflection  increases  as  for  large  K  and  is 
proportional  to  the  geometric  mean  of  the  individual  signal-to- 
noise  ratios. 

A  comparison  of  the  Gaussian  deflection  d^  in  (37)  with  some 
exact  results  for  the  desired  deflection  d_  in  (35)  will  be  made 
in  the  next  section  for  selected  values  of  K  and  A.  The 
asymptotic  behavior  (38)  will  not  be  employed  for  that 
comparison,  since  it  is  valid  only  for  larger  values  of  K. 
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EXACT  RESULTS 

Some  special  cases  for  mean  magnitude  ^J^  defined  in  (26)  can 
be  carried  out  in  closed  form.  These  results  complement  the 
earlier  approximations,  furnish  a  check  on  the  Gaussian 
approximations,  and  establish  their  regions  of  accuracy. 

MEAN  MAGNITUDE  FOR  ONE  PIECE,  K  =  1 

When  K  »  1,  the  magnitude  of  complex  estimate  G  follows 
immediately  from  (7)  as 

jj 

|G|  -  ((a^+aj)2  +  (p^+b^)2]  ((a^+c^)2  +  (p^+d^)^]  .  (39) 

The  desired  average  over  the  six  random  variables  involved  is 
conducted  by  first  holding  random  variables  and  Pj  fixed.  The 
conditional  average  over  the  remaining  four  random  variables  then 
factors  into  the  product  of  two  averages.  The  first  conditional 
average  can  be  expressed  as 
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J  O'’  ••  J 


de 


-n 


2no: 


exp 


r  - 


2  2 

2aj^r  COS0  -  2fi^r  sinG  + 


2o: 


r 

J  dr  exp 


_  £1  +  -^ 


2o" 


X 

(9), 

[2; 

^  *  * 

llj 


X  1*1 


F,(-.5;l;-R^)  ,  (40) 


(22).  But,  it  must  now  be  observed  that  no  random  variables 
remain  in  the  end  result,  meaning  that  no  further  averaging  is 
required!  A  similar  approach  can  be  used  for  the  second  term  in 
(39),  yielding  the  desired  exact  result  for  the  mean  of  |g|  as 


n 


Pj(K-l)  »  f  Oj,  Oy  iFj(-.5;l;-R^)  ^F^ (-.5; Ij-R^) 


(41) 


This  can  be  compared  with  the  corresponding  Gaussian 
approximation  according  to  (30)  and  (29),  namely 


»-lg(K.I)  .  1.**  Cy  (l+R^+Ry)**  1^1 


..  -..I.- 


(42) 


A  short  table  follows;  for  K  «  the  Gaussian  approximation  is 
anywhere  from  6%  to  13%  in  error,  over  this  range  of  values. 


Tabulation  of  Mean  Magnitude  for  K  i 


R 

R 

/^l(K=l) 

X 

y 

‘^x 

®x  "y 

Pl(K»l) 

0 

0 

1.57 

1.77 

1.13 

.5 

.25 

2.18 

2.43 

1.12 

.5 

.5 

2.40 

2.66 

1.11 

1 

.25 

2.55 

2.80 

1.10 

1 

.5 

2.81 

3.08 

1.10 

1 

1 

3.29 

3.56 

1.08 

2 

.5 

3.52 

3.77 

1.07 

2 

1 

4.12 

4.38 

1.06 
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DEFLECTION  FOR  ONE  PIECE,  K  -  1 

The  deflection  of  interest  was  defined  in  (35).  When  we  use 
exact  results  (41)  and  (34),  we  find,  that  for  K  **  1, 


dc(K-l) 


n 


(16  - 


^  [iri(-.5;l;-Rj^)  jFj(-.5;l;-Ry)  -  l] 


(43) 


On  the  other  hand,  the  Gaussian  approximation  to  the  deflection 
is  given  by  (37)  and  (29)  in  the  form 


<Jg(K-l) 


l*R,H.Ry. 


(44) 


A  comparison  of  these  two  results  is  given  in  the  following 
table;  the  Gaussian  approximation  overestimates  the  deflection  by 
about  40%  for  K  «  1.  This  is  not  too  surprising  when  we  recall 
that  the  Gaussian  approximation  cannot  be  expected  to  be  valid 
for  K  «  1,  but  rather  to  be  best  for  large  K,  where  summation 
variables  R  and  Q  in  (7)  are  tending  toward  Gaussian. 


Tabulation  of  Deflection  for  K  =  1 


^'x 

"y 

d^(K-l) 

dg(K=l) 

d^(K«l) 

0 

0 

0+ 

0+ 

1.51 

.5 

.25 

.489 

.707 

1.45 

.5 

.5 

.668 

.959 

1.43 

1 

.25 

.789 

1.11 

1.41 

1 

.5 

.999 

1.41 

1.41 

1 

1 

1.39 

1.93 

1.39 

2 

.5 

1.57 

2.16 

1.37 

2 

1 

2.06 

2.81 

1.37 
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JOINT  PROBABILITY  DENSITY  FUNCTION  FOR  SIGNAL  ABSENT,  ANY  K 


When  amplitude  A  in  input  signal  (1)  is  zero,  then  parameter 
Y  in  (9)  is  zero,  and  the  exact  joint  characteristic  function  of 
R  and  Q  in  (19)  reduces  to 

— K 

f°Q(i5,ih)  *  [l  +  8^(5^  +  r  (45) 


where  we  used  (15)  and  defined  B  «  a  o  /K.  The  corresponding 

X  y 

exact  joint  probability  density  function  of  R  and  Q  is  then 


Prq(U,V) 


dl  dY]  exp(-iu5-ivh) 


1  +  b" 


n 


0  -n 


de 


expr-ir(u  cos8  +  v  sinO)) 
(l  . 


Jo(pr) 


(l  +  B^r^) 


K 


K-1 

-B _ 

(K-1)! 


(46) 


2  2  k 

where  we  have  defined  p  *  (u  +  v  )  and  used  [2;  6.565  4].  The 

function  K^(x)  is  a  modified  Bessel  function  of  the  second  kind 

and  order  v  [3;  9.6].  Relation  (46),  which  applies  only  for 

signal  amplitude  A  »  0,  is  valid  for  all  argument  values  u,v  and 

any  number  of  pieces  K;  this  joint  density  is  seen  to  be  a 

2  2  k 

function  only  of  radius  (u  +  v  )^. 
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MOMENTS  OF  MAGNITUDE  ESTIMATE  |6|  FOR  SIGNAL  ABSENT,  ANY  K 

The  exact  mean  magnitude  >  |G|  is  given  by  (26).  When 
signal  amplitude  A  *  0,  this  becomes 


Pj(A-O)  -  JJ  du  dv  (u^  +  v^)**  p°q(u,v)  - 


J  dp  p  J  d+  p 


K-1 


-n 


n  2*“  (K-l)l  B 


K+1 


K  f^] 


T(K+.S) 
r(K 


(1/2) 

“x  “y  "  -jrrr  ““  *  ' 


(47) 


where  we  used  (46),  [2;  6.561  16],  and  B  ®  o  a  /K.  More 

X  y 

generally,  the  2v-th  moment  of  |Gi  for  signal  absent  is 
available  according  to 


JJ' 


du  dv  (u^+v^)’^  Prq(»'V) 


TfK-i-v)  r(v-t-l) 


r(K) 


2v 


(48) 


which  is  exact  for  all  K  and  v.  As  checks  on  (48),  we  have: 

2  2 

1  for  V  «  0;  result  (47)  for  v  »  %;  and  4  a  a  /K  for  v  »  1.  The 

X  y 

“2 

last  result  is  the  mean  square  value  R  +  Q  and  agrees  with 
exact  result  (34)  when  signal  amplitude  A  »  0  there.  The 
asymptotic  behavior  of  moment  (48)  is  given  by 


/i2v(A“0)  ~  r(v+l) 


4 


as  K  -»  «»  , 


(49) 


where  we  used  [3;  6.1.47]. 
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The  Gaussian  approximation,  Pjg(A*0),  to  is  given  by 

(31).  A  comparison  of  the  two  is  given  below;  it  reveals  that 
the  Gaussian  approximation  is  excellent  for  K  >  10. 


Tabulation  of  Mean  Magnitude  for  A  *  0 


K 

Pl(A~0)/(o^ay) 

Plg(A-0)/( 

1 

1.5708 

1.7725 

2 

1.1781 

1.2533 

3 

.9817 

1.0233 

4 

.8590 

.8862 

5 

.7731 

.7927 

6 

.7087 

.7236 

7 

.6581 

.6699 

8 

.6169 

.6267 

9 

.5827 

.5908 

10 

.5535 

.5605 

20 

.3939 

.3963 

30 

.3223 

.3236 

40 

.2794 

.2802 

50 

.2500 

.2507 

60 

.2283 

.2288 

70 

.2115 

.2118 

80 

.1979 

.1982 

90 

.1866 

.1868 

100 

.1770 

.1772 
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JOINT  PROBABILITY  DENSITY  FUNCTION  FOR  SIGNAL  PRESENT,  ANY  K 

We  have  analytically  determined  the  exact  joint 
characteristic  function  of  R  and  Q  in  (19),  for  arbitrary  signal 
amplitude  A  and  number  of  pieces  K.  However,  the  corresponding 
joint  probability  density  function  Pj^q(u,v)  is  not  generally 
available  in  closed  form.  When  (19)  is  substituted  into  the 
double  Fourier  transform  for  Pj^g(u,v),  and  a  change  to  polar 
coordinates  is  made,  the  following  single  integral  results: 


Prq(U,V) 


dr  r 


B(r) 


K 


exp 


2K  B(r] 

» 

jQ(D(r;u,v) ) 


(50) 


where 


®X  %  rf  V  21** 

B(r)  -  1  +  ^ -  ,  D(r;u,v)  -  r  [  [u  -  +  v*^]  .  (51) 

K 

Although  numerical  values  could  be  obtained  from  (50),  a  more 
efficient  approach  is  to  use  a  two-dimensional  fast  Fourier 
transform  directly  on  characteristic  function  (19). 

We  begin  by  defining,  for  numerical  convenience,  the 
normalized  random  variables 


The  joint  characteristic  function  of  r  and  q  is  then 


(52) 
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■  exp(i5r+ihq)  -  expf^p-^^ 
^  '  X  y 


RQlo^a„'<y^a„J 


X  y  X  Y 


r  2  2  21'^ 

1  +  (r+h  )/K^  exp  - 2LJt - 

^  '  1  +  (r 


U2rJrJ  -  (£^+h^)(VRy)/Kl 


+h^)/K^ 


.  (53) 


where  we  used  (19)  and  (22).  There  follows 

r  -  2  rJ  ,  q  -  0  ,  Tq  •  0  ,  oj  “  “  2(1+R  +R  )/K  .  (54) 

X  y  i  (j  X  y 

The  joint  probability  density  function  of  r  and  q  is 
p  (u,v)  -  JJ  d^  dh  exp(-iu^-ivh)  f  (i^^ih)  - 


m  « 


Re  J  d5  J  dh  exp(-iu^-ivti)  f  (i^/ih) 


— o  R®  2__.  2__. 

2x'^  k-0  X«-< 


exp[-ih^(uk+vX)]  fj.q(ikAj,iXZ\j]  ,  (55) 


where  we  used  the  conjugate  symmetry  of  took  a  common 

sampling  increment  for  f^g#  in  both  I  and  t).  Coefficient  Ej^ 
is  associated  with  the  trapezoidal  rule  and  is  1  for  all  k  except 
for  Eq  ■  1/2.  It  should  be  observed  that  the  resulting 
approximation  in  (55),  which  will  be  denoted  by  2j.g(u,v),  is 
periodic  in  both  u  and  v,  with  period  2n/^^,  This  aliased 
probability  density  function,  £j.q(u,v),  is  the  quantity  that  will 
be  evaluated. 
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We  now  take  samples  of  the  function  £j.q(u,v)  over  full 
periods  in  u  and  v,  that  is,  0  i.  m,n  i.  N-1,  according  to 


f2nm  2nn^ 
^rqlNAj'NAjJ 


2n 


00  00  /  •  O  \ 

Re  TIZ  exp[^=^^(mk+nX)  f  (ikA,,iXA 

k-O  X— »  ^  ;  rq  I 


aJ  N-1  N-1  .  . 

f(e  ZIZ  HZ  exppi^(mk+nX)  f  (ikA-,iXA,)  , 
k-O  X-O  ^  ;  a  I  r 


(56) 


where  {f^(ikA^,iXA^) )  is  the  collapsed  (or  prealiased)  version  of 
{Ek  f  j.q(ikA£,iXA£) } .  No  approximation  is  involved  in  the  last 
step  in  (56),  in  reducing  the  infinite  sums  to  finite  sums.  The 
double  sum  in  (56)  will  be  recognized  as  a  two-dimensional  fast 
Fourier  transform,  when  N  is  taken  as  a  power  of  2. 

The  common  increment  in  the  two  arguments  of  the  aliased 
probability  density  function  Ej.q(u,v)  in  (56)  is 


Sampling  increment  in  K  and  )i  must  be  small  enough  that  the 
resulting  aliasing  in  periodic  function  Ej.q(u,v)  is 
insignificant.  Also,  sampling  increment  A^  in  u  and  v  must  be 
small  enough  to  track  important  variations  in  This 

will  generally  require  large  values  of  N. 
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MEAN  MAGNITUDE  FOR  SIGNAL  PRESENT,  ANY  K 
The  desired  moment  is  mean  magnitude 

-  (R^  +  -  Ojj  Cy  (r^  +  - 

-  Oy  JJ  du  dv  (u^  +  Pj.q(u,v)  *  0^  Oy  D2  .  (58) 

In  order  to  evaluate  double  integral  D2,  three  approximations 
must  be  accepted.  First,  the  doubly  infinite  range  must  be 
replaced  by  a  square  of  size  2n/A^  covering  the  region  where 
p^q(u,v)  is  essentially  nonzero;  this  region,  to  be  denoted  by  S, 
is  roughly  centered  at  u,v  *  r,q.  Then,  Pj.q(U/V)  must  be 
replaced  by  Ej.q(u,v),  since  the  former  function  cannot  be 
evaluated.  Finally,  the  double  integral  must  be  replaced  by  a 
double  sum,  using  the  sample  points  furnished  by  (56).  The 
accuracy  of  these  three  replacements  depends  critically  on  the 
ability  to  accomplish  the  goals  listed  under  ( 57 ) ,  and  therefore 
on  the  ability  to  utilize  large  values  of  N  in  (56).  The 
resulting  approximation  to  D2  is 

D2  *  Ap  E  IZ  +  n^)**  Erq("4p>'>''p)  •  (5’l 

One  final  nuance  is  that  since  region  S  can  encompass 
negative  values  for  m  and/or  n,  whereas  (56)  is  typically 
evaluated  only  for  0  ^  m,n  N-1,  the  lookup  for  the  appropriate 
value  of  with  (m  +  n  )^  in  (59)  is  in  bins 

*  iBr  Mr 
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ID  modulo  N  and  n  modulo  N,  respectively.  A  program  that  achieves 
all  of  these  features  is  listed  in  the  appendix;  it  includes  some 
diagnostic  plots  that  keep  track  of  the  aliasing  and  attempt  to 
control  the  error  inherent  in  (59). 

An  example  for  K  -  1,  -  10  yielded  “  ^2  “ 

21.026487893,  when  done  exactly  by  means  of  (41).  As  an  illus¬ 
tration  of  the  accuracy  of  (59),  it  yielded  Dj  “  21.026487800  for 
the  same  parameter  values,  using  increment  «  .06  and  N  >  128. 
Also,  Kjij  ■  200  samples  of  characteristic  function  f^q(i^,i>))  in 
each  dimension  were  used,  thereby  minimizing  termination  error. 

Another  check  on  the  above  procedure  and  program  was 
accomplished  by  deliberately  taking,  as  a  test  case,  a  Gaussian 
two-dimensional  characteristic  function,  and  subjecting  it  to  the 
above  numerical  techniques.  The  exact  answer  for  the  mean 
magnitude  is  furnished  by  (30)  for  this  Gaussian  example.  In 
particular,  for  K  »  10,  R^  »  R^  ■  1,  (30)  yielded  A'ig/( Ojt’^y ^  “ 
2.1577687.  On  the  other  hand,  for  *  .6,  N  =  128,  «  50, 

numerical  procedure  (59)  yielded  2.1577675,  an  error  of  1.2E-6. 


DEFLECTION  OF  MAGNITUDE  |G| ,  ANY  K 

We  now  have  the  ability  to  exactly  evaluate  the  deflection  d 
of  magnitude  estimate  |G|  defined  in  (35),  and  to  compare  it  with 
the  Gaussian  approximation  defined  in  (36)  and  evaluated  in  (37). 
A  niunerical  comparison  is  presented  in  the  table  below. 
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Tabulation  of  Exact  and  Gaussian  Deflections 


K 

«x 

"y 

^1 

^ig 

^c 

9 

^f 

K 

m 

10 

.5 

.5 

1.2131 

1.2241 

2.1558 

2.2650 

.5 

50 

10 

1 

.5 

1.5978 

1.6068 

3.4135 

3.5712 

.5 

50 

10 

1 

1 

2.1515 

2.1578 

5.2231 

5.4517 

.5 

50 

10 

2 

.5 

2.1795 

2.1859 

5.3148 

5.5477 

.5 

50 

10 

2 

1 

2.9701 

2.9742 

7.8990 

8.2384 

.5 

50 

10 

2 

2 

4.1245 

4.1272 

11.672 

12.174 

.4 

50 

10 

4 

1 

4.1505 

4.1533 

11.757 

12.263 

.3 

50 

10 

4 

4 

8.1120 

8.1133 

24.706 

25.779 

.3 

50 

100 

1 

1 

2.0150 

2.0151 

19.748 

19.836 

1.5 

30 

For  K  >  10/  the  agreement  of  mean  magnitude  and  Gaussian 
approximation  is  very  good  over  the  entire  range  of  parameter 
values  considered,  with  the  Gaussian  approximation  being  a  slight 
overestimate  by  less  than  1%.  On  the  other  hand,  the  agreement 
between  deflections  d^  and  dg  is  not  quite  as  good,  with  the 
Gaussian  case  overestimating  by  about  4.5%.  High  accuracy  in  the 
deflection  for  K  »  10  can  only  be  achieved  through  the  detailed 
numerical  procedure  presented  above;  the  Gaussian  approximation 
has  some  limitations  at  this  low  value  of  K,  the  number  of 
independent  pieces. 

On  the  other  hand,  for  K  *  100,  the  means  are  virtually 
identical,  while  the  deflections  differ  by  0.5%.  This  is  an 
illustration  of  the  approach  of  summation  variables  R  and  Q  in 
(7)  to  Gaussian  for  large  numbers  of  pieces,  K. 
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AUTO  SPECTRUM  ESTIMATE 

In  this  section,  we  consider  adding  the  two  received 
processes  together  and  estimating  the  resulting  auto  spectrum. 

We  then  evaluate  the  deflection  of  this  auto  spectrum  estimate 
and  compare  it  with  the  deflection  for  the  magnitude  of  the 
complex  cross  spectrum  estimate,  |g| . 

Since  the  two  input  noises  x(t)  and  y(t)  utilized  in  (6)  and 
(7)  can  have  different  levels,  we  scale  them  and  sum  according  to 

2(t)  »  [s(t)  +  x(t)l  +  X[s(t)  +  y(t)]  = 

«  (1  +  X)  s(t)  +  x(t)  +  X  y(t)  .  (60) 

Scale  factor  X  will  be  chosen  to  maximize  the  deflection  of  the 
auto  spectrum  estimate  of  process  z(t).  (More  generally,  we 
should  filter  the  two  processes  and  add.) 


CHARACTERISTIC  FUNCTION  OF  AUTO  SPECTRUM  ESTIMATE 

Analogous  to  (4),  (5),  and  (6),  the  k-th  voltage  density 
estimate,  at  analysis  frequency  f,  of  process  z(t)  is 

Zj^  s  J  dt  exp(-i2iift)  Wj^(t)  z(t)  = 

=  (1  +  X)  I  +  X,^  +  X  - 

»  (1  +  X)(aj^  +  iPj^)  +  •  (61) 
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The  auto  spectrum  estimate  at  analysis  frequency  f  is  given  by 


K  K 

G,  IZki"  *  iOk 


)  +  e.  +  if^r  .  (62) 


where  independent  Gaussian  random  variables 


®k  “  *k  ^  *  ^k  ^  ‘^k  ' 

with  properties 

5;;  -  0  ,  1;;  -  0  ,  ej  -  .  (64) 

Here,  we  used  (12)  and  (13).  An  alternative  form  for  (62)  is 


((1  +  X)  + 


(65) 


We  now  hold  the  set  of  random  variables  {Oj^}  and  (Pj^)  fixed 
and  compute  the  conditional  characteristic  function  of  the  k-th 
term  of  (65).  Using  (64),  the  Gaussian  property  of  the  random 
variables  {ej^}  and  and  (17),  the  desired  quantity  is 
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But,  by  use  of  (9),  the  end  result  in  (66)  is  not  a  random 

variable  at  all,  and  furthermore,  does  not  depend  on  k. 

Therefore,  the  characteristic  function  of  auto  spectrum  estimate 

G,  in  (65)  is  given  by 
z 


T{U) 


exp 


(1  Y 

1  -  ij;2ai/K 

e 


(67) 


This  result  in  (67)  is  exact.  By  expanding  In  F(iC)  in  a 
power  series  in  i{;,  the  j-th  cumulant  of  estimate  is  found  to 
be 


(2.21’ 

(1  +  X)"  R  R  ] 

1  +  3  5 

X  R  +  R  J 

for  j  i  1 


(68) 


where  we  used  (64)  and  (22). 


DISTRIBUTION  OF  AUTO  SPECTRUM  ESTIMATE 


The  exceedance  distribution  function  corresponding  to 
characteristic  function  (67)  is  the  detection  probability  for 
random  variable  G_  and  is  given  by  [5] 

Z 


Pr(G,  >  V) 
z 


2K  (1  +  X)^  R^  Ry' 

K  V 

Js' 

X^  R  +  R 

X  y  ' 

9 

2  .  ,2  2 

la^  +  X  OyJ 

4 

(69) 


where  we  used  (64)  and  (22).  The  false  alarm  probability  is 

obtained  by  setting  R„  =  R„  *  0,  thereby  yielding 

X  y 
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P 


f 


exp 


Kv/2 


2  ^  ,2 
+  X 


(70) 


We  now  want  to  choose  scale  factor  X  in  (60)  so  as  to 
maximize  the  detection  probability  while  holding  the  false  alarm 
probability  fixed.  This  latter  requirement  means  holding  the 
argument  of  the  exponential  in  (70)  fixed,  which  makes  threshold 
V  a  function  of  X.  It  also  makes  the  second  argument  of  the 
function  in  (69)  constant.  Therefore,  maximization  of  is 
achieved  by  maximizing  the  first  argument  in  (69),  or 
equivalently  by  maximizing  the  quantity 


II  *  >^1^ 

+  R 


y 


by  choice  of  scale  factor  X.  The  best  choice  is 


(71) 


X 


(72) 


leading  to  maximum  value  1/R  +  1/R„  for  (71).  Substitution  of 

X  y 

these  results  in  (69)  yield  the  maximxim  detection  probability  as 


‘K 


h  o 

K  V 

_Z 

)  ' 

2  2 

a 

,a  +  a  , 

X 

xj 

(73) 


(The  nonsymmetry  in  the  second  argument  can  be  eliminated  by 
using  the  symmetric  combination  ^[s(t)  +  x(t)]/a^  +  Ji”^[s(t)  + 
y(t)]/ay  instead  of  (60),  and  choosing  I)  optimally.)  The 
corresponding  false  alarm  probability  is  obtained  by  replacing 
the  first  argument  in  (73)  by  zero,  that  is,  R  *  R  «  0. 

^  y 
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DEFLECTION  OF  AUTO  SPECTRUM  ESTIMATE 


The  deflection  of  auto  spectrum  estimate  6  is  defined 

Z 

analogously  to  that  for  magnitude  cross  spectrum  estimate  |6|  in 
(35),  namely 


Xi  -  Xi(A-O) 

X2{A-0)^ 


(1  +  X) 


R^  +  Ry 


(74) 


where  we  used  the  relevant  cumulants  in  (68).  But,  this  quantity 
has  exactly  the  same  dependence  on  scale  factor  X  as  does  (71). 
Therefore  the  best  choice  of  X  is  again  (72),  leading  to  the 
maximum  deflection,  which  is  exact  for  all  K,  of 

da  »  (Rjj  +  Ry)  .  (75) 


Thus,  the  choices  of  X  that  maximize  the  deflection  and  the 
detection  probability  coincide  for  the  auto  spectrum  estimate. 

The  maximiun  deflection  d^  in  (75)  is  not  always  larger  than 
the  magnitude  cross  spectrum  deflection  d  considered  earlier, 
even  though  (75)  has  utilized  the  best  scale  factor  X  in 
summation  (60).  For  example,  for  large  K,  we  have  from  (38)  the 
very  good  approximation, 

<*9  -  (4^  ■'x  "y)  as  K  -»  »  ■  (76 

The  ratio  of  deflections  is  therefore  given  by 


%  R. 


(77) 
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The  minimum  value  of  this  ratio  is  .926,  reached  when  R  *  R  . 

y  * 

However,  if  R„/R„  >  2.207  or  if  R„/R„  <  .453,  then  ratio  (77)  is 
y  X  y  X 

always  larger  than  1.  Thus,  which  deflection  is  larger  (for 

2  2 

large  K)  depends  on  the  ratio  R„/R„  *  small  K,  direct 

y  X  X  y 

numerical  evaluation  reveals  that  d.  is  usually  larger  than  d_. 

ot  C 


GRAPHICAL  COMPARISON  OF  DEFLECTIONS 

It  was  demonstrated  in  the  previous  section  that  the  Gaussian 
approximation  is  rather  accurate  for  evaluating  the  deflection  of 
the  magnitude  cross  spectrum  estimate  |g|,  when  K  is  larger  than 
10.  The  resulting  Gaussian  deflection  d^  was  given  by  (37)  along 
with  (29).  On  the  other  hand,  the  deflection  d^  for  the  auto 

a 

spectrum  estimate  6_  is  given  exactly  by  ( 75 ) ,  and  is  valid  for 
all  K. 

Plots  of  deflections  d^  and  d^^  are  presented  in  figures  1  and 

2  for  R„/R„  *  1  and  1/2,  respectively.  They  confirm  the  general 

y  * 

behavior  predicted  earlier.  For  example,  figure  1  for  R  *  R 

y  * 

shows  dg  to  be  larger  than  d^  for  large  K,  but  the  curves  cross 
for  smaller  values  of  K.  On  the  other  hand,  figure  2  for 
R  »  larger  than  d  ,  except  when  K  gets 

very  large.  The  ratio  R^/R^  =  1/2  is  not  smaller  than  the 

y  * 

breakpoint  .453  (above)  that  would  guarantee  d^  greater  than  dg 
for  large  K. 

The  fact  that  cross  spectrum  deflection  dg  (or  d^)  is  greater 
than  auto  spectznim  deflection  d^  for  some  ranges  of  the  parameter 
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TR  10709 


values  does  not  necessarily  reflect  in  the  relative  detection 
capability  of  the  two  processing  techniques.  After  all,  the 
deflection  criterion  only  involves  moments  up  through  second 
order,  whereas  the  full  detection  and  false  alarm  probabilities 
involve  all  orders  of  moments.  An  example  where  the  deflection 
£  a  random  variable  can  be  artificially  accentuated  is 
illustrated  in  the  next  section. 
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ACCENTUATION  OF  DEFLECTION 


The  deflection  of  a  random  variable  is  based  upon  its  two 
lowest  order  moments,  and  can  therefore  be  a  misleading  statistic 
regarding  detectability.  That  is,  the  detection  and  false  alarm 
probabilities  depend  on  the  entire  probability  density  functions 
for  signal  present  and  absent,  respectively,  not  just  their  first 
two  moments. 

To  illustrate  these  points,  consider  detection  of  a  Gaussian 
random  variable  z  with  mean  mQ  under  signal-absent  hypothesis  Hq, 
and  mean  m^^  (>  m^)  under  signal-present  hypothesis  H^.  Also,  let 
the  standard  deviations  have  a  common  value  a  under  both 
hypotheses.  Then,  the  deflection  of  random  variable  x  is 


“l  "  "*0 


(78) 


The  detection  probability,  for  threshold  v,  is  given  by 


-  Pr(x  >  v|Hj) 


00 

r  du 

(u  -  m^)^' 

1  k  ©Xp 

i  (2n)^  a 

o  2 

2a 

=  J  dt  (2n)”^  exp(-tV2)  s  *[-^)  •  (79) 

(v-mj^  )/a 

Similarly,  the  false  alarm  probability  is  given  by 

Pj  *  Pr(x  >  v|Hq)  =  •  (80) 

For  a  given  false  alarm  probability  P^,  (80)  can  be  solved  for 
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threshold  v,  and  then  substituted  into  (79).  The  result  is 

•■d  ■  ‘K  *  *<*■£))  '  <«1) 

where  f  is  the  inverse  i  function,  and  we  used  (78).  Thus,  given 
a  specified  performance  level  P^,  P^,  the  single  parameter  d^^ 
completely  quantifies  performance.  Observe  that  we  are  still 
using  the  entire  probability  density  functions  of  x  under  Hq  and 
H^,  as  we  must  in  order  to  evaluate  the  exceedance  distribution 
functions  in  (79)  and  (80);  however,  the  receiver  operating 
characteristic  depends  on  only  the  single  parameter  d^^,  through 
rule  (81). 

Now,  let  us  consider  a  monotonic  nonlinear  distortion  of 
random  variable  x,  yielding  new  random  variable  y  according  to 

y  -  exp(cf)  ,  (82) 

where  scaling  C  (>  0)  is  an  unspecified  constant  at  the  moment. 
Obviously,  the  receiver  operating  characteristic  for  random 
variable  y  will  be  identical  with  that  determined  for  x  above  in 
(79),  (80),  and  (81);  only  the  thresholds  will  change. 

However,  let  us  now  consider  the  deflection  of  random 
variable  y.  Since  x  is  Gaussian,  we  have  under  hypothesis  Hj^, 
the  n-th  moment  of  y  in  the  form 
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The  deflection  of  random  variable  y  follows  immediately  as 


exp(Cd^)  -  1 

dy  *  7 - 2 - * 

y  expCC*^)  -  1  ^ 


(84) 


Deflection  d  depends  on  d  and  the  dimensionless  scaling  C.  If 

y  ^ 

scaling  C  is  very  small,  then  we  have  d  ■  d  ;  this  agrees  with 

y  * 

the  observation  that  distortion  (82)  is  virtually  linear  then. 

However,  if  parameter  C  is  substantial,  deflection  d^  can  be 

much  greater  than  d^.  In  fact,  given  a  value  of  d^,  there  is  a 

value  of  C,  namely  C  =  d  ,  at  which  d  peaks,  with  value 

X  y 


max  d 

c  y 


(85) 


As  an  example,  the  value  of  d  is  greater  than  1000  if  d  >  3.72. 

y  * 

Thus,  the  deflection  of  random  variable  y  can  be  greatly 
accentuated  relative  to  the  deflection  of  x,  merely  by  performing 
a  monotonic  nonlinear  distortion.  The  ability  to  achieve  this 
artificial  improvement  in  deflection  strongly  cautions  against 
relying  on  the  deflection  as  a  reliable  measure  of  detectability. 
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SUMMARY 

The  joint  characteristic  function  of  the  real  and  imaginary 
parts  of  the  complex  cross  spectrum  estimate  has  been  derived  in 
closed  form,  for  arbitrary  signal  strength  and  noise  spectra. 

For  noise-only,  the  corresponding  joint  probability  density 
function  has  also  been  derived  in  closed  form  and  used  to  obtain 
exact  results  for  fractional  moments  of  the  magnitude  of  the 
cross  spectrum  estimate.  For  signal  present,  an  efficient  two- 
dimensional  fast  Fourier  transform  numerical  procedure  has  been 
utilized  to  get  accurate  probability  density  functions  and 
moments . 

When  the  number  of  pieces,  K,  used  in  the  estimate  of  the 
cross  spectrum  is  large,  a  Gaussian  approximation  has  been 
employed  for  the  joint  probability  density  function  of  the  real 
and  imaginary  parts.  Numerical  computations  reveal  that  this 
Gaussian  approximation  is  adequate  if  K  >  10,  and  is  very 
accurate  for  K  >  100.  This  Gaussian  approximation  has  then  been 
used  to  determine  the  deflection  of  the  magnitude  of  the  cross 
spectrum  estimate. 

Comparisons  of  the  deflections  for  the  magnitude  of  the  cross 
spectziun  estimate  and  for  the  auto  spectrum  estimate  reveal  that 
they  are  rather  close  to  each  other.  However,  even  though  one 
deflection  may  be  larger  than  the  other  for  some  ranges  of 
parameter  values,  that  does  not  necessarily  make  the 
corresponding  processor  a  better  detector.  An  example  is 
presented  to  show  how  the  deflection  may  be  artificially  enhanced 
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merely  by  nonlinear  transformation  of  the  decision  variable,  but 
without  any  change  in  the  fundamental  detectability  of  the 
signal. 

A  program  is  furnished  in  BASIC  which  enables  calculation  of 
the  joint  probability  density  function  of  the  real  and  imaginary 
parts  of  the  cross  spectrum  estimate,  for  arbitrary  signal 
strength.  In  addition,  it  calculates  the  mean  magnitude  of  the 
cross  spectrum  estimate  and  compares  it  with  the  Gaussian 
approx imat ion . 
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APPENDIX  -  PROGRAM  FOR  CALCULATION  OF  ^J^ 

This  appendix  contains  a  listing  of  a  BASIC  program  for  the 

evaluation  of  the  aliased  joint  probability  density  function 

given  by  (56),  in  addition  to  the  normalized  moment  D2  defined  by 

(59).  Inputs  required  of  the  user  are  Delf  in  line  10, 

P  (K)  in  line  20,  Rx  (R  )  in  line  30,  Ry  (R  )  in  line  40,  N  (N) 

X  y 

in  line  50,  and  Km  (K^)  in  line  60.  An  explanation  of  each  of 
these  symbols  is  given  in  the  program  listing. 

The  first  plot  produced  is  a  slice  of  the  magnitude  of  the 
aliased  characteristic  function  f_(i5,ifi)  in  (56)  along  the  I, 

OL 

axis;  this  affords  a  determination  of  whether  adequate  decay  has 
been  realized.  The  next  plot  displays  the  real  and  imaginary 
parts  of  fg^;  this  indicates  whether  the  sampling  rate  is  adequate 
to  track  the  variations  in  these  two  functions. 

Then,  the  sum  of  the  sampled  probability  density  function  is 
computed  and  subtracted  from  1;  this  error  furnished  a  measure  of 
the  accuracy  with  which  the  density  has  been  calculated.  Next,  a 
slice  of  density  2(u,v)  along  the  u  axis  is  plotted;  this 
indicates  whether  sufficient  decay  has  been  achieved  before  the 
aliasing  shows  up.  (Strictly,  this  observation  replaces  the  one 
above  on  the  real  and  imaginary  parts  of  )  The  user  must  also 
note  the  maximum  location  of  the  density  and  enter  this  number 
into  the  program  at  this  point.  Finally,  a  slice  of  p(u,v)  in  v, 
for  u  equal  to  the  maximum  location,  is  plotted;  this  guarantees 
that  adequate  decay  in  the  other  dimension  of  the  density  has 
been  achieved. 
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The  values  of  fj-^,  /Jj^(A=0),  d^,  and  dg  are  then  printed 

out.  This  complete  procedure  furnishes  a  measure  of  accuracy  of 
the  Gaussian  approximation  results,  provided  that  sufficient 
decays  have  been  realized  in  all  the  plots  indicated  above. 
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Delf*.5  !  INCREMENT  FOR  CHfiRfiCTER I  ST  I C  FUNCTION 

P=10.  !  NUMBER  OF  PIECES,  K 

Px=l.  !  MEfiSURE  OF  S I GNftL-TO-NO I SE  RATIO  IN  x 

Ry=l.  !  MEASURE  OF  SIGNAL-TO-NO I SE  RATIO  IN  y 

N=12S  !  SIZE  OF  FAST  FOURIER  TRANSFORM 

t<m=50  !  NUMBER  OF  SAMPLES  IN  EACH  DIMENSION 

DOUBLE  N,Km,Nl,K,L,Kl,Lt ,K1  !  INTEGERS 

DIM  Fr<127, 127>,Fi <127, 127),X<127),Y<127>,CosC32> 

N1=N-1 

REDIM  Fr<0;Nl,0:Nn,Fi<0;Nl,0;Nl),X<0;Nl),Y<0;Nl>,Cos<0:N/4> 
Tl*Delf*Delfx<P*P) 

T2®<Rx+Ry)*Del f ♦Del  f/P 
T3*2. ♦SQR<Rx^Ry)^De1 f 
R=2. *PI/N 
FOR  K=0  TO  N/'4 

Cos<K)=C0S<A^K)  !  QUARTER-COSINE  TABLE  IN  Cos<#) 

NEXT  K 

FOR  K=0  TO  Km 
Kt=K  MODULO  N 
K2*K*K 


T4=T3*K 

FOR  L*-Km  TO  Km 

Lt*L  MODULO  N 

Sq=K2+L*L 

T= 1 . +T 1 *Sq 

R=-P*L0G<T)-T2»Sq/T 

IF  A<-500.  THEN  320 

E=EXP<R) 

A=T4/T 

Fr <Kt , Lt )aFr<Kt , Lt )+E*C0S< A)  !  COLLAPSING 
Fi <Kt , Lt )=Fi <Kt , Lt )+E^SIN<R) 

NEXT  L 
NEXT  K 
GINIT 

GRAPHICS  ON 
UIND0I4  0,N,-10,0 
GRID  Nz-S,  1 
FOR  K=0  TO  N1 
Fr=Fr<K, 0) 

Fi aFi <K, 0) 

FsaFr*Fr+Fi ♦Fi 
IF  Fs>0.  THEN  450 
PENUP 
GOTO  460 

PLOT  K,LGT<Fs>^.5 

NEXT  K 

PENUP 

PRINT  "  I f <xi , 0)  I  ” 

PAUSE 

PRINT  “  Re  f<xi,0)  and  Im  fCxi.O)" 
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510  GCLERR 

520  MINDOU  0,N,-1, 1 

530  GRID  N/8, .2 

540  FOR  K«0  TO  N1 

550  PLOT  K,Fr<K,0)  !  Re  f(xi,0> 

560  NEXT  K 

570  PENUP 

580  LINE  TYPE  3 

590  FOR  Ks0  TO  Nl 

600  PLOT  K,Fi<K,0)  !  Im  f<xi,0) 

610  NEXT  K 

620  PENUP 

630  LINE  TYPE  1 

640  FOR  K«0  TO  Nl 

650  FOR  L»0  TO  Nl 

660  X<L)»Fr<K,L) 

670  Y<L)*Fi<K,L) 

680  NEXT  L 

690  IF  K>0  THEN  720 

700  MAT  X*X*<.5) 

710  MAT  Y»Y*<.5) 

720  CALL  Fftl4<N,Cos<e),X<*: ,Y<*>) 

730  .  FOR  L=0  TO  Nl 

740  Fr<K,L)=X<L) 

750  Fi<K,L)=Y<L) 

760  NEXT  L 

770  NEXT  K 

780  FOR  L=0  TO  Nl 

790  FOR  K=0  TO  Nl 

800  X<K)=Fr<K,L) 

810  Y<K)«Fi<K,L) 

820  NEXT  K 

830  CALL  Fftl4<N,Cos<*),X<»: ,Y<»)) 

840  FOR  K=0  TO  Nl 

850  Fr<K,L)=X<K) 

860  !  Fi<K,L)=Y<K)  !  Fi <*)  unnecessary 

870  NEXT  K 

880  NEXT  L 

890  MAT  Fr=Fr*<Delf*Delf/'<2.*PI*PI)) 

900  Del  p=2.  ♦PI/'CNeDel  f  ) 

910  S=SUM<Fr)*Delp*Delp 

920  PRINT  "P  =";P;''  Rx  =';Rx;"  Ry  =";Ry 

930  PRINT  "ERROR  *";S-1. 

940  PRINT  “  p<u,0)“ 

950  GCLEAR 

960  UINDOM  0,N,-12,O 

970  GRID  N/8,1 

980  FOR  K=0  TO  Nl 

990  Fr»Fr<K,0) 

1000  IF  Fr<>0.  THEN  1030 

1010  PENUP 

1020  GOTO  1040 

1030  PLOT  K,LGT<RBS<Fr)) 

1040  NEXT  K 

1050  PLOT  N,LGT<Fr<0,0)) 

1060  PENUP 

1070  INPUT  "MAXIMUM  LOCATION  OF  Fr<K,0>;",Kl 

1080  PRINT  "  p<ui,u)’ 

1090  GCLEAR 

1100  UINDOU  0,N,-12,e 
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1 110  GRID  N/8,  1 

1120  FOR  L=0  TO  N1 

1130  Fr=Fr<Kl,L) 

1140  IF  Fr<>0.  THEN  1170 

1150  PENUP 

1160  GOTO  1180 

1170  PLOT  L,LGT<flBS<Fr)) 

1180  NEXT  L 

1190  PLOT  N.LGTCFrCKl.O)) 

1200  PENUP 

1210  PRUSE 

1220  GCLERR 

1230  ril»0. 

1240  FOR  K*Kl-N/2  TO  Kl+N/2 

1250  Kt-K  MODULO  N 

1260  K2=K*K 

1270  FOR  L»-N/'2  TO  N-^2 

1280  Lt*L  MODULO  N 

1290  Ml=Ml+SQR<t^2+L*L)*Fr<Kl ,  Lt  ) 

1300  NEXT  L 

1310  NEXT  K 

1320  Ml=Ml*Del p*De1 p*De1 p 

1330  PRINT  Mul  =";M1 

1340  V=P*Rx»Ry.^<  1 . +Rx  +  Ry> 

1350  CRLL  Fll<-.5, 1. ,-V,Fll,I > 

1360  Mlg=SQR<Pl»<l.+Rx+Ry)/P)*Fll 

1370  PRINT  “MuUGRUSS)  =''5Mlg 

1380  M10sPI 

1390  FOR  K=1  TO  P 

1400  M10=M10*<K-.5)/K 

1410  NEXT  K 

1420  PRINT  "Mul<0)  =";M10 

1430  M2  =  4.  'P 

1440  Dc  =  <Ml-M10>/SQR<M2-M10*rM0) 

1450  PRINT  "  Dc  ='';Dc 

1460  Dg  =  SQR<PI/<4.-Pn)*<SQR(l.+Rx+Ry'»*Fll-l.  ) 

1470  PRINT  ■•Dc<GflUSS)  =";Dg 

1480  PRINT 

1490  PRUSE 

1500  END 

1510  ! 

152v'  SUB  Fftl4<D0UBLE  N.RERL  Cos<*)  ,  X<*)  ,  Y<*)  )  !  N<  =2^^  1 4=  1 6384 ;  0  SUBS 
1530  DOUBLE  Log2n ,  N 1 ,  N2 ,  N3 ,  N4  ,  J ,  K  !  INTEGERS  <  2-^31  =  2,147,483,648 
1540  DOUBLE  II, 12, 13, 14, 15, 16, 17, 18, 19, 110, Ill, 112, 113, 114, LCOilS) 

1550  IF  N=1  THEN  SUBEXIT 

1560  IF  N>2  THEN  1640 

1570  R=XC0)+X<1) 

1580  X( 1 )=X<0)-X< n 

1590  X(0)=R 

1600  R»Y<0)+Y<1) 

1610  Y<1)»Y<0>-Ya> 

1620  Y<0)=R 

1630  SUBEXIT 

1640  R=L0G<N)/L0G<2. ) 

1650  Log2n-R 
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1660 

1670 

1680 

1690 

1700 

1710 

1720 

1730 

1740 

1750 

1760 

1770 

1780 

1790 

1800 

1810 

1820 

1830 

1840 

1850 

1860 

1870 

1880 

1890 

1900 

1910 

1920 

1930 

1940 

1950 

I960 

1970 

1980 

1990 

2000 

2010 

2020 

2030 

2040 

2050 

2060 

2070 

2080 

2090 

2100 

2110 

2120 

2130 

2140 

2150 

2160 

2170 

2180 

2190 

2200 


IF  flBS<fi-Log2n)<l.E-8  THEN  1690 

PRINT  "N  »”;N;"IS  NOT  fl  POWER  OF  2;  DISALLOWED." 

PAUSE 

Nl»N/4 

N2=N1+1 

N3-N2+1 

N4-N3+N1 

FOR  Il»l  TO  Log2n 
I2»2''<Log2n-Il) 

13*2*12 

I4«N/'I3 

FOR  15*1  TO  12 

I6*<I5-1)*I4+1 

IF  I6<*N2  THEN  1830 

Al*-Cos<N4-I6-U 

A2*-Cos<I6-Nl-l) 

GOTO  1850 
Al*Cos<I6-l) 

A2*-Cos<N3-I6-l) 

FOR  17*0  TO  N-I3  STEP  13 

18*17+15-1 

19*18+12 

T1*X<I8) 

T2*X< 19) 

T3=Y< 18) 

T4*Y<I9) 

A3*T1-T2 

A4ST3-T4 

X<I8)«T1+T2 

Y<I8)*T3+T4 

X<I9)=A1*A3-A2*A4 

Y<I9)*A1*A4+A2*A3 

NEXT  17 

NEXT  15 

NEXT  II 

I l*Log2n+l 

FOR  12*1  TO  14 

L<I2-1)*1 

IF  I2>Log2n  THEN  2060 
L<I2-l)s2''<Il-I2) 

NEXT  12 
K*0 

FOR  11*1  TO  L<13) 

FOR  12*11  TO  L<12)  STEP  L<13) 

FOR  13*12  TO  L<11)  STEP  L<12) 

FOR  14*13  TO  L<10)  STEP  L<11) 

FOR  15*14  TO  L<9)  STEP  L<10) 

FOR  16*15  TO  L<8)  STEP  L<9) 

FOR  17*16  TO  L<7)  STEP  L<8) 

FOR  18*17  TO  L<6)  STEP  L<7) 

FOR  19*18  TO  L<5)  STEP  L<6) 

FOR  110*19  TO  L<4)  STEP  L<5) 

FOR  111*110  TO  L<3)  STEP  L<4) 

FOR  112*111  TO  L<2)  STEP  LCS) 

FOR  113*112  TO  LCD  STEP  L<2) 
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2210 

2220 

2230 

2240 

2250 

2260 

2270 

2280 

2290 

2300 

2310 

2320 

2330 

2340 

2350 

2360 

2370 

2380 

2390 

2400 

2410 

2420 

2430 

2440 

2450 

2460 

2470 

2480 

2490 

2500 

2510 

2520 

2530 

2540 

2550 

2560 

2570 

2580 

2590 

2600 

2610 

2620 

2630 

2640 

2650 

2660 

2670 

2680 

2690 

2700 

2710 

2720 

2730 

2740 

2750 

2760 


FOR  I14»I13  TO  L<0)  STEP  L<1) 

J-I14-1 

IF  K>J  THEN  2300 
fl*X<K) 

XCK)»X<J) 

X<J)»fl 

fi»Y<K) 

Y<K)«Y<J) 

Y<J>«fl 
K«K+1 
NEXT  114 
NEXT  113 
NEXT  112 
NEXT  Ill 
NEXT  110 
NEXT  19 
NEXT  18 
NEXT  17 
NEXT  16 
NEXT  15 
NEXT  14 
NEXT  13 
NEXT  12 
NEXT  II 
SUBEND 
! 

SUB  Fll<fi,B,X,Fll,D)  !  POWER  SERIES 

Error=l.E-16  !  RELATIVE  ERROR  TOLERANCE 

Numbersl000  !  MAXIMUM  NUMBER  OF  TERMS  IN  SERIES 

DOUBLE  Number, N  !  INTEGERS 

B1»B-1. 

IF  X<0.  THEN  2640 
A1=A-1. 

FI l»T»Bi g=l . 

FOR  N«1  TO  Number 
Fn*FLT<N) 

T*T*X*<Fn+Al)/'<Fn*<Fn+Bl)) 

F11=F11+T 

Af»ABS<Fll) 

Big=MAX<Big,Rf ) 

IF  ABS<T)<=Error*Af  THEN  2750 
NEXT  N 
GOTO  2740 
Bal»B-A-l. 

Fll=T=Big=EXP<X) 

FOR  N=1  TO  Number 
Fn»FLT<N) 

T=-T*X*<Fn+Bal)/'<Fn*<Fn+Bl)) 

F11»F11+T 

Af«ABS<Fll) 

B1g«MAX(Big,Af > 

IF  ABS<:T)<=Error*Af  THEN  2750 
NEXT  N 

PRINT  Number; "TERMS  IN  SUB  Fll  AT  ";A;B;X 
D»15.'LGT<Big/Af )  ‘  NUMBER  OF  SIGNIFICANT  DIGITS 

SUBEND 
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